不走寻常路的单细胞表达量矩阵

最近看到了一个单细胞数据集，是GSE133283，标题是：《Single cell transcriptome analysis of non-neuronal cell populations in mouse brain cortex》，蛮有意思的。在主页：https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE133283

可以看到11个样品都是10X单细胞技术的， 并且作者给出来了文件：

GSM3904816_Adult-1_gene_counts.tsv.gz 32.0 Mb
GSM3904817_Adult-2_gene_counts.tsv.gz 32.2 Mb
GSM3904818_Adult-3_gene_counts.tsv.gz 32.7 Mb
GSM3904819_Aged-1_gene_counts.tsv.gz 37.0 Mb
GSM3904820_Aged-2_gene_counts.tsv.gz 37.0 Mb
GSM3904821_EAE-1_gene_counts.tsv.gz 64.1 Mb
GSM3904822_EAE-2_gene_counts.tsv.gz 63.4 Mb
GSM3904823_Juvenile-1_gene_counts.tsv.gz 29.5 Mb
GSM3904824_Juvenile-2_gene_counts.tsv.gz 29.3 Mb
GSM3904825_Juvenile-3_gene_counts.tsv.gz 44.2 Mb
GSM3904826_Juvenile-4_gene_counts.tsv.gz 40.0 Mb 

我本来以为是表达量矩阵，准备使用如下所示的代码：

###### step1:导入数据 ######   
library(data.table)
dir='GSE135045_RAW/' 
samples=list.files( dir )
samples
# samples = head(samples,10)

library(data.table)
sceList = lapply(samples,function(pro){ 
  # pro=samples[1] 
  print(pro) 
  ct=fread( file.path(dir,pro),data.table = F)
  ct[1:4,1:4]
  ct[nrow(ct),1:4]
  rownames(ct)=ct[,1]
  ct=ct[,-1]
  ct[1:4,1:4]
  sce =CreateSeuratObject(counts =  ct , 
                          min.cells = 5,
                          min.features = 300 )
  return(sce)
})
names(sceList)  

但是发现居然每个txt文件读入，都在   ct[1:4,1:4] 就失败了，我肉眼看了看，居然只有3 列！！！查看其中一个样品，如下所示：

>   head(ct)
           gene             cell count
1 0610005C13Rik AACCGCGTCCGTTGCT     1
2 0610005C13Rik ATTATCCAGTTGTCGT     1
3 0610005C13Rik CAGTAACAGCAGCCTC     1
4 0610005C13Rik CATCGAAGTAGCGTAG     1
5 0610005C13Rik CCACTACTCACATGCA     1
6 0610005C13Rik CCTAAAGCAAGTAATG     1

如果大家单细胞经验丰富，就会发现它其实有点类似于 sparse Matrix of class "dgCMatrix" ，只不过呢，它删除了全部的基因表达量为0的情况，如果2万个基因和2千个细胞本来是一个有2千万个值，但是单细胞里面90%~95%都是0 ，所以其实就只剩下了  20~200万个值，所以 ，我们简单的把数据形式 长变成宽即可：

>  dim(ct)
[1] 6182813       3
>   library(reshape2)
>   tmp = dcast(ct,gene~cell,value.var = 'count')
>   dim(tmp)
[1] 21892  3664
>   ct[1:4,1:4]
Error in `[.data.frame`(ct, 1:4, 1:4) : 选择了未定义的列
>   ct = tmp
>   ct[1:4,1:4]
           gene AAACCTGAGATGTGTA AAACCTGAGGTACTCT AAACCTGAGTGTTAGA
1 0610005C13Rik               NA               NA               NA
2 0610007N19Rik               NA               NA               NA
3 0610007P14Rik                1                2               NA
4 0610009B14Rik               NA               NA               NA

接下来就是同样的简单的代码：

###### step1:导入数据 ######   
library(data.table)
dir='GSE133283_RAW' 
samples=list.files( dir )
samples
# samples = head(samples,10)
library(data.table)
sceList = lapply(samples,function(pro){ 
  # pro=samples[1] 
  print(pro) 
  ct=fread( file.path(dir,pro),data.table = F)
  head(ct)
  library(reshape2)
  tmp = dcast(ct,gene~cell,value.var = 'count')
  dim(tmp)
  ct = tmp
  ct[1:4,1:4]
  ct[nrow(ct),1:4]
  rownames(ct)=ct[,1]
  ct=ct[,-1]
  ct[1:4,1:4]
  ct[is.na(ct)]=0
  ct[1:4,1:4]
  sce =CreateSeuratObject(counts =  ct ,
                          project =   gsub('_gene_counts.tsv.gz','',gsub('^GSM[0-9]*_','',pro) ) ,
                          min.cells = 5,
                          min.features = 300 )
  return(sce)
})
names(sceList)  

最后就是走标准的降维聚类分群即可：

标准的降维聚类分群

需要一些背景知识给这些单细胞亚群进行命名：

脑部单细胞亚群各自基因

往期回顾

---

---

---

如果你对单细胞转录组研究感兴趣，但又不知道如何入门，也许你可以关注一下下面的课程

看完记得顺手点个“在看”哦！

生物 | 单细胞 | 转录组丨资料每天都精彩

长按扫码可关注